using FreeBird
walkers = AtomWalker.(generate_initial_configs(120, 562.5, 6))
lj = LJParameters(epsilon=0.1, sigma=2.5)LJParameters(0.1 eV, 2.5 Å, Inf, 0.0 eV)Engineering a computational surface modeling toolbox in Julia
Ray M. Yang
Washington University in St. Louis
2024-11-04
@ PChem Group Meeting
Overview
Introduction to
FreeBird.jl
Design, functionality, and extensibilityIntroduction to Research Software Engineering (RSE)
What is RSE, and how does it fit in Chemistry
What is FreeBird.jl
- It’s a computational package
- It’s developed by the Wexler Group
- It does surface free energy calculations via nested sampling
- It’s written in Julia, a modern, high-level programming language
What is Julia?
is a high-level, high-performance, dynamic programming language for technical computing. It’s fast, easy to learn, and has a syntax that is familiar to users of other technical computing environments.
- Version 1.0 was released in 2018
- Developed by mathematician Alan Edelman et al. at MIT
- It’s designed for high-performance numerical analysis and computational science
- It’s free and open source
FreeBird.jl design
- It’s modular
- It’s extensible
- It’s user-friendly
- It’s well-documented
- It’s open-source
FreeBird.jl functionality
- Systems: atomistic and lattice models
- Energy calculator: Lennard-Jones and lattice Hamiltonians
- Methods: nested sampling, Metropolis Monte Carlo, exact enumeration, Wang-Landau sampling
Nested sampling
Atomistic walkers
Atomistic walkers
Atomistic walkers
FastSystem(H₆, periodic = FFF):
bounding_box : [ 15 0 0;
0 15 0;
0 0 15]u"Å"
AtomView(H, [ 4.01502, 4.18315, 10.6918]u"Å")
AtomView(H, [ 4.94146, 3.03619, 4.64721]u"Å")
AtomView(H, [ 5.35381, 11.0872, 7.85456]u"Å")
AtomView(H, [ 9.05429, 1.99427, 7.42124]u"Å")
AtomView(H, [ 13.382, 13.0819, 0.368069]u"Å")
AtomView(H, [ 13.3007, 11.1216, 0.182911]u"Å")
.------------------------------------.
/| |
/ | |
/ | |
/ | |
/ | |
/ | |
/ | |
/ | |
* | |
| | |
| | H |
| | H |
| | |
| | |
| | |
| | |
| | H |
| .------------------------------------.
| / H /
| / H H /
| / /
| / /
| / /
| / /
| / /
|/ /
*------------------------------------* Construct the liveset
LJAtomWalkers(AtomWalker{1}, LJParameters):
[1] AtomWalker{1}(
configuration : FastSystem(H₆, periodic = FFF, bounding_box = [[14.999999999999998, 0.0, 0.0], [0.0, 14.999999999999998, 0.0], [0.0, 0.0, 14.999999999999998]]u"Å")
energy : 5.2747930136178764 eV
iter : 0
list_num_par : [6]
frozen : Bool[0]
energy_frozen_part : 0.0 eV)
[2] AtomWalker{1}(
configuration : FastSystem(H₆, periodic = FFF, bounding_box = [[14.999999999999998, 0.0, 0.0], [0.0, 14.999999999999998, 0.0], [0.0, 0.0, 14.999999999999998]]u"Å")
energy : -0.03877740224016672 eV
iter : 0
list_num_par : [6]
frozen : Bool[0]
energy_frozen_part : 0.0 eV)
[3] AtomWalker{1}(
configuration : FastSystem(H₆, periodic = FFF, bounding_box = [[14.999999999999998, 0.0, 0.0], [0.0, 14.999999999999998, 0.0], [0.0, 0.0, 14.999999999999998]]u"Å")
energy : -0.04244190870170766 eV
iter : 0
list_num_par : [6]
frozen : Bool[0]
energy_frozen_part : 0.0 eV)
[4] AtomWalker{1}(
configuration : FastSystem(H₆, periodic = FFF, bounding_box = [[14.999999999999998, 0.0, 0.0], [0.0, 14.999999999999998, 0.0], [0.0, 0.0, 14.999999999999998]]u"Å")
energy : -0.10913780161856197 eV
iter : 0
list_num_par : [6]
frozen : Bool[0]
energy_frozen_part : 0.0 eV)
[5] AtomWalker{1}(
configuration : FastSystem(H₆, periodic = FFF, bounding_box = [[14.999999999999998, 0.0, 0.0], [0.0, 14.999999999999998, 0.0], [0.0, 0.0, 14.999999999999998]]u"Å")
energy : -0.11327569913084067 eV
iter : 0
list_num_par : [6]
frozen : Bool[0]
energy_frozen_part : 0.0 eV)
⋮
Omitted 110 walkers
⋮
[116] AtomWalker{1}(
configuration : FastSystem(H₆, periodic = FFF, bounding_box = [[14.999999999999998, 0.0, 0.0], [0.0, 14.999999999999998, 0.0], [0.0, 0.0, 14.999999999999998]]u"Å")
energy : -0.07630973295538208 eV
iter : 0
list_num_par : [6]
frozen : Bool[0]
energy_frozen_part : 0.0 eV)
[117] AtomWalker{1}(
configuration : FastSystem(H₆, periodic = FFF, bounding_box = [[14.999999999999998, 0.0, 0.0], [0.0, 14.999999999999998, 0.0], [0.0, 0.0, 14.999999999999998]]u"Å")
energy : 1.3074885636307922 eV
iter : 0
list_num_par : [6]
frozen : Bool[0]
energy_frozen_part : 0.0 eV)
[118] AtomWalker{1}(
configuration : FastSystem(H₆, periodic = FFF, bounding_box = [[14.999999999999998, 0.0, 0.0], [0.0, 14.999999999999998, 0.0], [0.0, 0.0, 14.999999999999998]]u"Å")
energy : 132.93929690640402 eV
iter : 0
list_num_par : [6]
frozen : Bool[0]
energy_frozen_part : 0.0 eV)
[119] AtomWalker{1}(
configuration : FastSystem(H₆, periodic = FFF, bounding_box = [[14.999999999999998, 0.0, 0.0], [0.0, 14.999999999999998, 0.0], [0.0, 0.0, 14.999999999999998]]u"Å")
energy : 31.242106865621327 eV
iter : 0
list_num_par : [6]
frozen : Bool[0]
energy_frozen_part : 0.0 eV)
[120] AtomWalker{1}(
configuration : FastSystem(H₆, periodic = FFF, bounding_box = [[14.999999999999998, 0.0, 0.0], [0.0, 14.999999999999998, 0.0], [0.0, 0.0, 14.999999999999998]]u"Å")
energy : -0.009495483209446563 eV
iter : 0
list_num_par : [6]
frozen : Bool[0]
energy_frozen_part : 0.0 eV)
LJParameters(0.1 eV, 2.5 Å, Inf, 0.0 eV)Run nested sampling
Nested sampling results
FastSystem(H₆, periodic = FFF):
bounding_box : [ 15 0 0;
0 15 0;
0 0 15]u"Å"
AtomView(H, [ 2.92958, 6.58633, 3.60833]u"Å")
AtomView(H, [ 2.1724, 6.20473, 0.94135]u"Å")
AtomView(H, [0.227731, 6.37336, 2.94226]u"Å")
AtomView(H, [ 3.21941, 8.66147, 1.76325]u"Å")
AtomView(H, [0.510936, 8.44946, 1.10138]u"Å")
AtomView(H, [ 1.269, 8.83089, 3.7619]u"Å")
.------------------------------------.
/| |
/ | |
/ | |
/ | |
/ | |
/ | |
/ | |
/ | |
* | |
| | |
| | |
| | |
| | |
| | |
| | |
| | |
| H| |
| .------------------------------------.
| / H /
| H / H /
| H /
| / H /
| / /
| / /
| / /
|/ /
*------------------------------------* Nested sampling results
Nested sampling results
Run nested sampling
Nested sampling in 7 lines of code!
ls = LJAtomWalkers(AtomWalker.(generate_initial_configs(120, 562.5, 6)), LJParameters())
mc = MCRandomWalkClone()
ns_params = NestedSamplingParameters(200, 0.1, 0.01, 1e-5, 1.0, 0, 200)
save = SaveEveryN(n_traj=10, n_snap=10_000)
energies, liveset, _ = nested_sampling_loop!(ls, ns_params, 10_000, mc, save)Or, in a single line if you really want to!
Accessing the documentation
help?> nested_sampling_loop!
search: nested_sampling_loop! nested_sampling_step!
nested_sampling_loop!(liveset::AtomWalkers, ns_params::NestedSamplingParameters, n_steps::Int64, mc_routine::MCRoutine; args...)
Perform a nested sampling loop for a given number of steps.
Arguments
≡≡≡≡≡≡≡≡≡
• liveset::AtomWalkers: The initial set of walkers.
• ns_params::NestedSamplingParameters: The parameters for nested sampling.
• n_steps::Int64: The number of steps to perform.
• mc_routine::MCRoutine: The Monte Carlo routine to use.
Returns
≡≡≡≡≡≡≡
• df: A DataFrame containing the iteration number and maximum energy for each step.
• liveset: The updated set of walkers.
• ns_params: The updated nested sampling parameters.
───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
nested_sampling_loop!(liveset::LatticeGasWalkers, ns_params::LatticeNestedSamplingParameters, n_steps::Int64, mc_routine::MCRoutine; args...)
Perform a nested sampling loop on a lattice gas system for a given number of steps.
Arguments
≡≡≡≡≡≡≡≡≡
• liveset::LatticeGasWalkers: The initial set of walkers.
• ns_params::LatticeNestedSamplingParameters: The parameters for nested sampling.
• n_steps::Int64: The number of steps to perform.
• mc_routine::MCRoutine: The Monte Carlo routine to use.
Keyword Arguments
≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡
• args...: Additional arguments.
Returns
≡≡≡≡≡≡≡
• df: A DataFrame containing the iteration number and maximum energy for each step.
• liveset: The updated set of walkers.
• ns_params: The updated nested sampling parameters.
───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
nested_sampling_loop!(liveset::AtomWalkers, n_steps::Int64, mc_routine::MixedMCRoutine, save_strategy::DataSavingStrategy)
Perform a nested sampling loop for a given number of steps.
Arguments
≡≡≡≡≡≡≡≡≡
• liveset::AtomWalkers: The initial set of walkers.
• n_steps::Int64: The number of steps to perform.
• mc_routine::MixedMCRoutine: The mixed Monte Carlo routine to use.
• save_strategy::DataSavingStrategy: The strategy for saving data.
Returns
≡≡≡≡≡≡≡
• df::DataFrame: The data frame containing the iteration number and maximum energy for each step.
• liveset::AtomWalkers: The updated set of walkers.
• mc_routine.ns_params_main: The updated nested sampling parameters for the main routine.Accessing the documentation
help?> NestedSamplingParameters
search: NestedSamplingParameters LatticeNestedSamplingParameters nested_sampling_step!
mutable struct NestedSamplingParameters <: SamplingParameters
The NestedSamplingParameters struct represents the parameters used in the nested sampling scheme.
Fields
≡≡≡≡≡≡
• mc_steps::Int64: The number of total Monte Carlo moves to perform.
• initial_step_size::Float64: The initial step size, which is the fallback step size if MC routine fails to accept a move.
• step_size::Float64: The on-the-fly step size used in the sampling process.
• step_size_lo::Float64: The lower bound of the step size.
• step_size_up::Float64: The upper bound of the step size.
• fail_count::Int64: The number of failed MC moves in a row.
• allowed_fail_count::Int64: The maximum number of failed MC moves allowed before resetting the step size.Finding phase transitions
# Calculate the ω factors
ωi = ωᵢ(energies.iter, 120)
# Shift the energies to be greater than or equal to zero
Ei = energies.emax .- minimum(energies.emax)
# Specify the temperatures that we are interested in
Ts = collect(1:0.1:1000)
# Define the Boltzmann constant
kb = 8.617333262e-5 # eV/K
# Calculate the inverse temperatures
β = 1 ./(kb.*Ts)
# Define the degrees of freedom, which is 3×6 for the 6-particle system
dof = 18
# Calculate the partition functions for each temperature
Zs = [partition_function(b, ωi, Ei) for b in β]
# Calculate the internal energies for each temperature
U = [internal_energy(b, ωi, Ei) for b in β]
# Calculate the heat capacities as a function of temperature
cvs = cv(energies, β, dof, 120)9991-element Vector{Float64}:
0.0013175732622637696
0.0013163656401052968
0.0013164865839072522
0.0013180887545059053
0.0013209878572451813
0.0013247504441241784
0.0013288500299832545
0.0013327866085318337
0.0013361503619112659
0.0013386436487862145
0.0013400792149547592
0.001340367194218547
0.001339497919186227
⋮
0.000865730605363293
0.0008657001748975948
0.0008656697600629232
0.0008656393608486451
0.0008656089772439815
0.0008655786092383255
0.0008655482568209475
0.000865517919981241
0.0008654875987085314
0.0008654572929921685
0.0008654270028215192
0.0008653967281859931Finding surface phase transitions
Finding surface phase transitions
Surface stuff
using AtomsBase
configs = read_configs("slab.extxyz",pbc="TTF")
configs = FastSystem.(configs)
slabs = AtomWalker{2}.(configs, list_num_par=[80,6],frozen=[1,0])480-element Vector{AtomWalker{2}}:
AtomWalker{2}(
configuration : FastSystem(H₈₆, periodic = TTF, bounding_box = [[11.224620483093732, 0.0, 0.0], [0.0, 9.72080648619833, 0.0], [0.0, 0.0, 29.164864246657352]]u"Å")
energy : 0.0 eV
iter : 0
list_num_par : [80, 6]
frozen : Bool[1, 0]
energy_frozen_part : 0.0 eV)
AtomWalker{2}(
configuration : FastSystem(H₈₆, periodic = TTF, bounding_box = [[11.224620483093732, 0.0, 0.0], [0.0, 9.72080648619833, 0.0], [0.0, 0.0, 29.164864246657352]]u"Å")
energy : 0.0 eV
iter : 0
list_num_par : [80, 6]
frozen : Bool[1, 0]
energy_frozen_part : 0.0 eV)
AtomWalker{2}(
configuration : FastSystem(H₈₆, periodic = TTF, bounding_box = [[11.224620483093732, 0.0, 0.0], [0.0, 9.72080648619833, 0.0], [0.0, 0.0, 29.164864246657352]]u"Å")
energy : 0.0 eV
iter : 0
list_num_par : [80, 6]
frozen : Bool[1, 0]
energy_frozen_part : 0.0 eV)
AtomWalker{2}(
configuration : FastSystem(H₈₆, periodic = TTF, bounding_box = [[11.224620483093732, 0.0, 0.0], [0.0, 9.72080648619833, 0.0], [0.0, 0.0, 29.164864246657352]]u"Å")
energy : 0.0 eV
iter : 0
list_num_par : [80, 6]
frozen : Bool[1, 0]
energy_frozen_part : 0.0 eV)
AtomWalker{2}(
configuration : FastSystem(H₈₆, periodic = TTF, bounding_box = [[11.224620483093732, 0.0, 0.0], [0.0, 9.72080648619833, 0.0], [0.0, 0.0, 29.164864246657352]]u"Å")
energy : 0.0 eV
iter : 0
list_num_par : [80, 6]
frozen : Bool[1, 0]
energy_frozen_part : 0.0 eV)
AtomWalker{2}(
configuration : FastSystem(H₈₆, periodic = TTF, bounding_box = [[11.224620483093732, 0.0, 0.0], [0.0, 9.72080648619833, 0.0], [0.0, 0.0, 29.164864246657352]]u"Å")
energy : 0.0 eV
iter : 0
list_num_par : [80, 6]
frozen : Bool[1, 0]
energy_frozen_part : 0.0 eV)
AtomWalker{2}(
configuration : FastSystem(H₈₆, periodic = TTF, bounding_box = [[11.224620483093732, 0.0, 0.0], [0.0, 9.72080648619833, 0.0], [0.0, 0.0, 29.164864246657352]]u"Å")
energy : 0.0 eV
iter : 0
list_num_par : [80, 6]
frozen : Bool[1, 0]
energy_frozen_part : 0.0 eV)
AtomWalker{2}(
configuration : FastSystem(H₈₆, periodic = TTF, bounding_box = [[11.224620483093732, 0.0, 0.0], [0.0, 9.72080648619833, 0.0], [0.0, 0.0, 29.164864246657352]]u"Å")
energy : 0.0 eV
iter : 0
list_num_par : [80, 6]
frozen : Bool[1, 0]
energy_frozen_part : 0.0 eV)
AtomWalker{2}(
configuration : FastSystem(H₈₆, periodic = TTF, bounding_box = [[11.224620483093732, 0.0, 0.0], [0.0, 9.72080648619833, 0.0], [0.0, 0.0, 29.164864246657352]]u"Å")
energy : 0.0 eV
iter : 0
list_num_par : [80, 6]
frozen : Bool[1, 0]
energy_frozen_part : 0.0 eV)
AtomWalker{2}(
configuration : FastSystem(H₈₆, periodic = TTF, bounding_box = [[11.224620483093732, 0.0, 0.0], [0.0, 9.72080648619833, 0.0], [0.0, 0.0, 29.164864246657352]]u"Å")
energy : 0.0 eV
iter : 0
list_num_par : [80, 6]
frozen : Bool[1, 0]
energy_frozen_part : 0.0 eV)
AtomWalker{2}(
configuration : FastSystem(H₈₆, periodic = TTF, bounding_box = [[11.224620483093732, 0.0, 0.0], [0.0, 9.72080648619833, 0.0], [0.0, 0.0, 29.164864246657352]]u"Å")
energy : 0.0 eV
iter : 0
list_num_par : [80, 6]
frozen : Bool[1, 0]
energy_frozen_part : 0.0 eV)
AtomWalker{2}(
configuration : FastSystem(H₈₆, periodic = TTF, bounding_box = [[11.224620483093732, 0.0, 0.0], [0.0, 9.72080648619833, 0.0], [0.0, 0.0, 29.164864246657352]]u"Å")
energy : 0.0 eV
iter : 0
list_num_par : [80, 6]
frozen : Bool[1, 0]
energy_frozen_part : 0.0 eV)
AtomWalker{2}(
configuration : FastSystem(H₈₆, periodic = TTF, bounding_box = [[11.224620483093732, 0.0, 0.0], [0.0, 9.72080648619833, 0.0], [0.0, 0.0, 29.164864246657352]]u"Å")
energy : 0.0 eV
iter : 0
list_num_par : [80, 6]
frozen : Bool[1, 0]
energy_frozen_part : 0.0 eV)
⋮
AtomWalker{2}(
configuration : FastSystem(H₈₆, periodic = TTF, bounding_box = [[11.224620483093732, 0.0, 0.0], [0.0, 9.72080648619833, 0.0], [0.0, 0.0, 29.164864246657352]]u"Å")
energy : 0.0 eV
iter : 0
list_num_par : [80, 6]
frozen : Bool[1, 0]
energy_frozen_part : 0.0 eV)
AtomWalker{2}(
configuration : FastSystem(H₈₆, periodic = TTF, bounding_box = [[11.224620483093732, 0.0, 0.0], [0.0, 9.72080648619833, 0.0], [0.0, 0.0, 29.164864246657352]]u"Å")
energy : 0.0 eV
iter : 0
list_num_par : [80, 6]
frozen : Bool[1, 0]
energy_frozen_part : 0.0 eV)
AtomWalker{2}(
configuration : FastSystem(H₈₆, periodic = TTF, bounding_box = [[11.224620483093732, 0.0, 0.0], [0.0, 9.72080648619833, 0.0], [0.0, 0.0, 29.164864246657352]]u"Å")
energy : 0.0 eV
iter : 0
list_num_par : [80, 6]
frozen : Bool[1, 0]
energy_frozen_part : 0.0 eV)
AtomWalker{2}(
configuration : FastSystem(H₈₆, periodic = TTF, bounding_box = [[11.224620483093732, 0.0, 0.0], [0.0, 9.72080648619833, 0.0], [0.0, 0.0, 29.164864246657352]]u"Å")
energy : 0.0 eV
iter : 0
list_num_par : [80, 6]
frozen : Bool[1, 0]
energy_frozen_part : 0.0 eV)
AtomWalker{2}(
configuration : FastSystem(H₈₆, periodic = TTF, bounding_box = [[11.224620483093732, 0.0, 0.0], [0.0, 9.72080648619833, 0.0], [0.0, 0.0, 29.164864246657352]]u"Å")
energy : 0.0 eV
iter : 0
list_num_par : [80, 6]
frozen : Bool[1, 0]
energy_frozen_part : 0.0 eV)
AtomWalker{2}(
configuration : FastSystem(H₈₆, periodic = TTF, bounding_box = [[11.224620483093732, 0.0, 0.0], [0.0, 9.72080648619833, 0.0], [0.0, 0.0, 29.164864246657352]]u"Å")
energy : 0.0 eV
iter : 0
list_num_par : [80, 6]
frozen : Bool[1, 0]
energy_frozen_part : 0.0 eV)
AtomWalker{2}(
configuration : FastSystem(H₈₆, periodic = TTF, bounding_box = [[11.224620483093732, 0.0, 0.0], [0.0, 9.72080648619833, 0.0], [0.0, 0.0, 29.164864246657352]]u"Å")
energy : 0.0 eV
iter : 0
list_num_par : [80, 6]
frozen : Bool[1, 0]
energy_frozen_part : 0.0 eV)
AtomWalker{2}(
configuration : FastSystem(H₈₆, periodic = TTF, bounding_box = [[11.224620483093732, 0.0, 0.0], [0.0, 9.72080648619833, 0.0], [0.0, 0.0, 29.164864246657352]]u"Å")
energy : 0.0 eV
iter : 0
list_num_par : [80, 6]
frozen : Bool[1, 0]
energy_frozen_part : 0.0 eV)
AtomWalker{2}(
configuration : FastSystem(H₈₆, periodic = TTF, bounding_box = [[11.224620483093732, 0.0, 0.0], [0.0, 9.72080648619833, 0.0], [0.0, 0.0, 29.164864246657352]]u"Å")
energy : 0.0 eV
iter : 0
list_num_par : [80, 6]
frozen : Bool[1, 0]
energy_frozen_part : 0.0 eV)
AtomWalker{2}(
configuration : FastSystem(H₈₆, periodic = TTF, bounding_box = [[11.224620483093732, 0.0, 0.0], [0.0, 9.72080648619833, 0.0], [0.0, 0.0, 29.164864246657352]]u"Å")
energy : 0.0 eV
iter : 0
list_num_par : [80, 6]
frozen : Bool[1, 0]
energy_frozen_part : 0.0 eV)
AtomWalker{2}(
configuration : FastSystem(H₈₆, periodic = TTF, bounding_box = [[11.224620483093732, 0.0, 0.0], [0.0, 9.72080648619833, 0.0], [0.0, 0.0, 29.164864246657352]]u"Å")
energy : 0.0 eV
iter : 0
list_num_par : [80, 6]
frozen : Bool[1, 0]
energy_frozen_part : 0.0 eV)
AtomWalker{2}(
configuration : FastSystem(H₈₆, periodic = TTF, bounding_box = [[11.224620483093732, 0.0, 0.0], [0.0, 9.72080648619833, 0.0], [0.0, 0.0, 29.164864246657352]]u"Å")
energy : 0.0 eV
iter : 0
list_num_par : [80, 6]
frozen : Bool[1, 0]
energy_frozen_part : 0.0 eV)Surface stuff
FastSystem(H₈₆, periodic = TTF):
bounding_box : [ 11.2246 0 0;
0 9.72081 0;
0 0 29.1649]u"Å"
.---------------------------.
/| |
/ | |
/ | |
/ | |
/ | |
* | |
| | |
| | |
| | |
| | |
| | |
| | |
| | |
| | |
| | |
| | H H |
| | |
| | |
| | |
| | H H |
| | H |
| | |
| | |
| | |
| | |
| | H H H H |
| | H |
| H | HH HH HH H |
| H| H H H |
| H | H H H |
H HHH HHH HHH HH |
| | |
|H |HH HH HH H |
| H | H H H |
| | H H H H |
| HH .--HHH----HHH----HHH----H---.
| / /
| HH HH HH HH /
H / HH HH HH H /
| / /
|H H H H /
*---------------------------* Surface stuff
surf_lj = LJParameters(epsilon=0.1, sigma=2.5, cutoff=4.0)
surf_ls = LJAtomWalkers(slabs, surf_lj)
mc = MCRandomWalkClone()
surf_ns_params = NestedSamplingParameters(500, 0.1, 0.01, 1e-5, 1.0, 0, 200)
save = SaveEveryN(n_traj=10, n_snap=50_000)
surf_energies, liveset, _ = nested_sampling_loop!(surf_ls, surf_ns_params, 50_000, mc, save)(50000×2 DataFrame Row │ iter emax │ Int64 Float64 ───────┼───────────────────── 1 │ 1 5.71085e12 2 │ 2 4.59391e9 3 │ 3 4.44588e9 4 │ 4 7.43776e8 5 │ 5 4.21399e8 6 │ 6 8.31856e7 7 │ 7 2.29957e7 8 │ 8 2.0287e7 9 │ 9 1.1382e7 10 │ 10 8.46159e6 11 │ 11 7.13711e6 ⋮ │ ⋮ ⋮ 49991 │ 49991 -58.5274 49992 │ 49992 -58.5274 49993 │ 49993 -58.5274 49994 │ 49994 -58.5274 49995 │ 49995 -58.5274 49996 │ 49996 -58.5274 49997 │ 49997 -58.5274 49998 │ 49998 -58.5274 49999 │ 49999 -58.5274 50000 │ 50000 -58.5274 49979 rows omitted, LJAtomWalkers(AtomWalker{2}, LJParameters): [1] AtomWalker{2}( configuration : FastSystem(H₈₆, periodic = TTF, bounding_box = [[11.224620483093732, 0.0, 0.0], [0.0, 9.72080648619833, 0.0], [0.0, 0.0, 29.164864246657352]]u"Å") energy : -58.52736668197494 eV iter : 50000 list_num_par : [80, 6] frozen : Bool[1, 0] energy_frozen_part : -54.68859657133728 eV) [2] AtomWalker{2}( configuration : FastSystem(H₈₆, periodic = TTF, bounding_box = [[11.224620483093732, 0.0, 0.0], [0.0, 9.72080648619833, 0.0], [0.0, 0.0, 29.164864246657352]]u"Å") energy : -58.527366976199524 eV iter : 50000 list_num_par : [80, 6] frozen : Bool[1, 0] energy_frozen_part : -54.68859657133728 eV) [3] AtomWalker{2}( configuration : FastSystem(H₈₆, periodic = TTF, bounding_box = [[11.224620483093732, 0.0, 0.0], [0.0, 9.72080648619833, 0.0], [0.0, 0.0, 29.164864246657352]]u"Å") energy : -58.52736700308327 eV iter : 50000 list_num_par : [80, 6] frozen : Bool[1, 0] energy_frozen_part : -54.68859657133728 eV) [4] AtomWalker{2}( configuration : FastSystem(H₈₆, periodic = TTF, bounding_box = [[11.224620483093732, 0.0, 0.0], [0.0, 9.72080648619833, 0.0], [0.0, 0.0, 29.164864246657352]]u"Å") energy : -58.52736700465829 eV iter : 50000 list_num_par : [80, 6] frozen : Bool[1, 0] energy_frozen_part : -54.68859657133728 eV) [5] AtomWalker{2}( configuration : FastSystem(H₈₆, periodic = TTF, bounding_box = [[11.224620483093732, 0.0, 0.0], [0.0, 9.72080648619833, 0.0], [0.0, 0.0, 29.164864246657352]]u"Å") energy : -58.527367183158844 eV iter : 50000 list_num_par : [80, 6] frozen : Bool[1, 0] energy_frozen_part : -54.68859657133728 eV) ⋮ Omitted 470 walkers ⋮ [476] AtomWalker{2}( configuration : FastSystem(H₈₆, periodic = TTF, bounding_box = [[11.224620483093732, 0.0, 0.0], [0.0, 9.72080648619833, 0.0], [0.0, 0.0, 29.164864246657352]]u"Å") energy : -58.527535068052806 eV iter : 50000 list_num_par : [80, 6] frozen : Bool[1, 0] energy_frozen_part : -54.68859657133728 eV) [477] AtomWalker{2}( configuration : FastSystem(H₈₆, periodic = TTF, bounding_box = [[11.224620483093732, 0.0, 0.0], [0.0, 9.72080648619833, 0.0], [0.0, 0.0, 29.164864246657352]]u"Å") energy : -58.52753654030334 eV iter : 50000 list_num_par : [80, 6] frozen : Bool[1, 0] energy_frozen_part : -54.68859657133728 eV) [478] AtomWalker{2}( configuration : FastSystem(H₈₆, periodic = TTF, bounding_box = [[11.224620483093732, 0.0, 0.0], [0.0, 9.72080648619833, 0.0], [0.0, 0.0, 29.164864246657352]]u"Å") energy : -58.52753791580169 eV iter : 50000 list_num_par : [80, 6] frozen : Bool[1, 0] energy_frozen_part : -54.68859657133728 eV) [479] AtomWalker{2}( configuration : FastSystem(H₈₆, periodic = TTF, bounding_box = [[11.224620483093732, 0.0, 0.0], [0.0, 9.72080648619833, 0.0], [0.0, 0.0, 29.164864246657352]]u"Å") energy : -58.527562374940956 eV iter : 50000 list_num_par : [80, 6] frozen : Bool[1, 0] energy_frozen_part : -54.68859657133728 eV) [480] AtomWalker{2}( configuration : FastSystem(H₈₆, periodic = TTF, bounding_box = [[11.224620483093732, 0.0, 0.0], [0.0, 9.72080648619833, 0.0], [0.0, 0.0, 29.164864246657352]]u"Å") energy : -58.52744314188147 eV iter : 50000 list_num_par : [80, 6] frozen : Bool[1, 0] energy_frozen_part : -54.68859657133728 eV) LJParameters(0.1 eV, 2.5 Å, 4.0, -9.763240814208984e-5 eV) , NestedSamplingParameters(500, 0.1, 0.009679022355186586, 1.0e-5, 1.0, 0, 200))
Surface results
FastSystem(H₈₆, periodic = TTF):
bounding_box : [ 11.2246 0 0;
0 9.72081 0;
0 0 29.1649]u"Å"
.---------------------------.
/| |
/ | |
/ | |
/ | |
/ | |
* | |
| | |
| | |
| | |
| | |
| | |
| | |
| | |
| | |
| | |
| | |
| | |
| | |
| | |
| | |
| | |
| | |
| H| |
| |H H |
| | |
| | H H HH HH |
| | H |
| H | HH HH HH H |
| H| H H H |
| H | H H H |
H HHH HHH HHH HH |
| | |
|H |HH HH HH H |
| H | H H H |
| | H H H H |
| HH .--HHH----HHH----HHH----H---.
| / /
| HH HH HH HH /
H / HH HH HH H /
| / /
|H H H H /
*---------------------------* Surface results
Surface results
# Calculate the ω factors
ωi = ωᵢ(surf_energies.iter, 480)
# Shift the energies to be greater than or equal to zero
Ei = surf_energies.emax .- minimum(surf_energies.emax)
# Specify the temperatures that we are interested in
Ts = collect(1:0.1:1000)
# Define the Boltzmann constant
kb = 8.617333262e-5 # eV/K
# Calculate the inverse temperatures
β = 1 ./(kb.*Ts)
# Define the degrees of freedom, which is 3×6 for the 6-particle system
dof = 18
# Calculate the partition functions for each temperature
Zs = [partition_function(b, ωi, Ei) for b in β]
# Calculate the heat capacities as a function of temperature
surf_cvs = cv(surf_energies, β, dof, 480)9991-element Vector{Float64}:
0.0013842345024584256
0.0014037564760710597
0.001422094212288646
0.0014393991765428245
0.0014549553127525684
0.0014679912279158577
0.0014781013099614134
0.0014853322720637684
0.0014900807662374433
0.0014929308416672571
0.0014945083183555742
0.001495382151241112
0.0014960136907975394
⋮
0.0030352937093422747
0.003035254740728882
0.0030352155242118813
0.003035176059862044
0.003035136347750496
0.003035096387948459
0.0030350561805273834
0.003035015725558813
0.0030349750231146062
0.0030349340732666397
0.003034892876087085
0.003034851431648211Surface results
Lattice systems
using Distributions, DataFrames
square_supercell_dimensions = (4, 4, 1)
square_basis=[(0.0, 0.0, 0.0)]
adsorption_energy = -0.04
nn_energy = -0.01
nnn_energy = -0.0025
# Initialize the lattice
occupied_sites = sample(1:16, 6, replace=false)
initial_lattice = LatticeSystem{SquareLattice}(;
supercell_dimensions = square_supercell_dimensions,
occupations=occupied_sites)
h = GenericLatticeHamiltonian(adsorption_energy, [nn_energy, nnn_energy], u"eV")GenericLatticeHamiltonian{2,Quantity{Float64, 𝐋² 𝐌 𝐓⁻², Unitful.FreeUnits{(eV,), 𝐋² 𝐌 𝐓⁻², nothing}}}:
on_site_interaction: -0.04 eV
nth_neighbor_interactions: [-0.01, -0.0025] eVLattice systems
walkers = [deepcopy(initial_lattice) for i in 1:2000]
for walker in walkers
walker.occupations = [false for i in 1:length(walker.occupations)]
for i in sample(1:length(walker.occupations), 6, replace=false)
walker.occupations[i] = true
end
end
unique!(x -> x.occupations, walkers) # remove duplicates
# Now, we can create the liveset by combining the walkers and the potential
ls = LatticeGasWalkers(deepcopy(LatticeWalker.(walkers[1:1000])), h, perturb_energy=1e-10)
# We can now create a Monte Carlo object
mc = MCNewSample()
# Let's specify the nested sampling parameters:
# Monte Carlo steps, energy perturbation, current fail count, allowed fail count
# Use `?LatticeNestedSamplingParameters` to see the documentation for more information
ns_params = LatticeNestedSamplingParameters(1,1e-10,0,500_000)
# We can also create a save object, using default value for most parameters except `n_snap` for saving snapshots
save = SaveEveryN(n_traj=10000, n_snap=500_000)
# And finally, we can run the simulation
lattice_energies, ls, _ = nested_sampling_loop!(ls, ns_params, 10_000, mc, save; accept_same_config=true)Lattice systems
ωi = ωᵢ(lattice_energies.iter, 1000)
# Shift the energies to be greater than or equal to zero
Ei = lattice_energies.emax .- minimum(lattice_energies.emax)
# Specify the temperatures that we are interested in
Ts = collect(1:0.1:200)
# Define the Boltzmann constant
kb = 8.617333262e-5 # eV/K
# Calculate the inverse temperatures
β = 1 ./(kb.*Ts)
# Define the degrees of freedom, which is 2×4 for the 4-particle system
dof = 8
# Calculate the partition functions for each temperature
Zs = [partition_function(b, ωi, Ei) for b in β]
# Calculate the heat capacities as a function of temperature
lattice_cvs = cv(lattice_energies, β, dof, 1000)1991-element Vector{Float64}:
0.00034469333056180083
0.00034469333142480505
0.0003446933376293437
0.00034469336959996606
0.0003446934965557751
0.0003446939063884268
0.0003446950258828272
0.0003446976939299924
0.0003447033750094083
0.00034471438632263944
0.00034473410463458824
0.00034476711998688886
0.0003448193116499381
⋮
0.00036374954549554323
0.0003637334434403411
0.0003637173611604454
0.0003637012986247458
0.00036368525580218975
0.0003636692326617828
0.0003636532291725878
0.00036363724530372545
0.00036362128102437375
0.00036360533630376795
0.00036358941111120074
0.0003635735054160216Lattice results
What’s next
- Publish! There are some journals that are specifically designed for software papers
Or, traditional journals like:
Take-home message
- Version control is a powerful tool
- It’s a good practice for code development
- It’s essential for collaborative projects
- It has great importance in scientific reproducibility
Start using it, now!
